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Magnetism in 2D atomic sheets has attracted considerable interest as its existence 
could allow the development of electronic and spintronic devices. The existence of 
magnetism is not sufficient for devices, however, as states must be addressable and 
modifiable through the application of an external drive. We show that defects in 
hexagonal boron nitride present a strong interplay between the the N-N distance 
in the edge and the magnetic moments of the defects. By stress-induced geometry 
modifications, we change the ground state magnetic moment of the defects. This 
control is made possible by the triangular shape of the defects as well as the strong 
spin localisation in the magnetic state. 
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Studied for only a few years, two-dimensional atomic sheets have shown a remarkably 
rich set of properties that could play a significant role in the development of spintronics and 
of the next generation of electronic devices. In particular, hexagonal boron nitride (h-BN) 
sheets have gained significant interest recently due to their similarities with graphene 1 ' 2 , 
both being 2D honeycomb lattices. However, h-BN differs from the latter on a number 
of properties. For example, it presents a wide bandgap associated with the partially ionic 
character of its bonding in contrast with the semi-metallic character of graphene 3 . As BN 
materials are generally more resistant to oxidation than graphene, h-BN should also be more 
suitable for applications at high-temperatures 4 . Its binary composition also increases the 
possibility for creating and positioning defects with specific properties, opening the door to a 
very pointed design and a rapid incorporation into technological devices. Recent studies, for 
example, have revealed interesting features of defects created by electron beam irradiation of 
h-BN sheets 5 " 8 . High-resolution transmission electron microscopy (HRTEM) shows triangu- 
lar multivacancy structures with discrete sizes but with a unique orientation. By means of 
the exit-wave reconstruction method in combination with HRTEM images, Jin et a/. 5 were 
able to assign a chemical label to the atoms surrounding the defects and demonstrate that 
these triangular multivacancies have N-terminated zigzag edges. This feature suggests that 
the boron atoms are easier to knock out than nitrogen under electron irradiation. This bias 
could facilitate the manufacturing of functional devices with well-controlled defect struc- 
tures and fully justifies a more careful study of these defects' properties. Du et al 9 have 
showed that these triangular vacancies exhibit a total magnetic moment proportional to the 
number of nitrogen dangling bonds in the N-terminated edges. This result is in line with 
calculations that predict large magnetic moments for BN nanoribbons with a ferromagnetic 
orientation for the N-zigzag edge and an antiferromagnetic for the B-edge 10 . While these 
numerical results have not yet received an experimental confirmation, they are supported 
by a characterization of defect-induced paramagnetism in graphene that demonstrates the 
validity of similar calculations in 2D carbon 11 . It is therefore tempting to see applications 
of h-BN in the current electronic industry but also for spintronic devices. For this purpose, 
however, static magnetic moments are not really sufficient, and tunable spin states remain 
the target. 

In this letter, we focus on the magnetic properties of multivacancy holes in h-BN sheets 
and their relation to structural features of the holes. Using spin-polarised density func- 
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tional (DFT) calculations as implemente in BigDFT 12 , we have characterized in detail the 
electronic state of the small multivacancy defects and explored approaches to tune their 
spin state. We have focused more particularly on the destabilisation of the reconstructed 
corner under tensile strain. Our results show that, under a given uniaxial deformation, 
the atoms surrounding a specific type of triangular multivacancy defect can undergo a pro- 
nounced geometric reconstruction. This reconstruction strongly affects the local electronic 
states, modifying the defect's total magnetic moment. These results reveal a tight interplay 
between the local strain and the local spin state, as the electronic states remain strongly 
localised during the transformation. This opens new ways to exploit the magnetic properties 
of defects in BN sheets. 

Calculations were done in a cell size that guarantees negligible elastic interactions between 
the periodic defects. To this end we use an orthorhombic supercell (35.2 x 34.8 A 2 ) with 
448 atoms and only the T point. In the perpendicular direction we use free boundary 
conditions, i.e. infinite vacuum. The wavelet grid is chosen such that the energy per atom is 
accurate within 5 meV, which corresponds to a grid spacing of 0.35 bohrs and a localisation 
radius of 9.6 bohrs. Geometries are considered optimized when the forces on atoms are 
less than 15 me V/A. We use a PBE exchange-correlation functional together with HGH 
pseudopotentials 13 . 

As starting point we address the formation energy of multivacancies in h-BN. This calcu- 
lation requires the use of a suitable chemical potential Hn(hb) for N(B) atoms. Because of 
the wide bonding diversity for both elements, we select, as point of comparison, compounds 
with an sp 2 environment similar to the h-BN sheet. Thus, knowing the chemical potential 
He in graphene, fijsr is obtained with the help of a sheet of g-C3N 4 14 , which is a crystalline 
graphitic carbon nitride compound. Hence, the formation energy for a vacancy defect can 
be estimated from the expression: 

E Dt = E x - D ? - (E x - n B ii B - n N fi N ), (1) 

where E x ~ Dt is the total energy of the h-BN sheet with a defect D T and E x is the total 
energy of the perfect system. In the following we consider several triangular multivacancies 
defects Vr, where Ub and number of B and N atoms have been removed, such as T = 
\Jtlb + tin. The correponding defects of size T 2 are either with N (V?) or with B (Vj?) 
edges. Calculated formation energies for T = 1, 2, 3 are displayed in Fig. 1(a). is found 
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(a) T T 2 (size) E(v") E(V^) 



1 1 6.71 8.17 

2 4 12.05 17.16 

3 9 19.93 23.96 




FIG. 1. (Color online) (a) Formation energies (in eV) for different N-edge and B-edge terminated 
vacancy sizes T 2 (T = 1, 2, 3). The values are given per defect, (b) Growth route for V^, V<± and 
N-terminated bigger holes. B atoms are in red and N atoms in blue. 

to be more favourable than V 8 . However, the incident electron energy in the experiments 5 ' 6 
is rather large with respect to the emission energy threshold for both atoms 15 . Thus, it is 
likely that both V-^ and V® are formed during the early irradiation time, while the former is 
the dominant one. Our results explain this discrepancy and offer a growth route for bigger 
defects. The direct formation of the defect (12.05 eV) is unfavoured with respect to 
an indirect formation starting from V 8 (+3.88 eV), see Fig. 1(b). This would explain the 
observed low concentration of Vj 8 in favor of V<^ . Conversely, the formation of the V<^ defect 
is unlikely from both routes, direct («~ 17 eV) and indirect, Vf* — > 10 eV). This could 

explain the large presence of in experimental observations rather than Vf . The next 
step, the formation of V^, is done by removing one biline (Fig. 1(b)) from one edge of 
for a cost of only 7.9 eV. The further transformation, ^(t+i)> 1S then expected to 

occur with an even lower energy cost. This proposed growth route is in agreement with the 
main experimental observations: a low concentration of V 8 \ a high concentration of Vf* and 
N-terminated big holes. The validity of this hypothesis could be easily assessed by varying 
the irradiation time and following the evolution of the defects. Of particular importance is 
the fact that by tuning the irradiation time it would be possible to restrict the defect family 
to Vi*, V 8 and V^, while long irradiation time always lead to large holes 5 ' 16 . 

The calculated total magnetic moment for V^,V^ and defects is 1, and 3 
respectively, and doesn't follow Lieb's theorem 17 contrary to graphene 18 . This violation is 
mainly due to strong reconstructions. Indeed, calculations show that the Vjf defects, with 
T > 1, undergo a significant reconstruction with the formation of a N-N bond in each corner 
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FIG. 2. (Color online), (a) Optimized geometry and Wannier centers of a h-BN sheet with one 
defect. B atoms are in red, N atoms are in blue and Wannier centers are in green. Zigzag 
(ZZ) and armchair (AC) directions are shown. The Wannier states outside the triangle are not 
disturbed by the defect, (b) Total energy and N-N distance at each corner (inset) as a function 
of the imposed magnetic moment, (c) Optimized geometry for the system with M — 2 \±b along 
with the difference in the spin magnetic charge density (p^ — pjj projected on the plane of the 
monolayer. For clearness, only atoms along the borders of the vacancy are included. 

of the triangle. Figure 2(a) shows the ground state geometry of a defect along with 
the Wannier centers from a Wannier projection 19 of the electronic states. The N-atoms at 
the corners form pairs with an average distance of only 1.66 A, 0.85 A shorter than in the 
perfect h-BN sheet. This pairing satisfies all dangling bonds, producing a homopolar a bond 
between the two N atoms which can be seen by the presence of only one Wannier center in 
the middle of the pair. This bonding reduces the total energy as well as the total magnetic 
moment associated with the hole to zero. The electronic states modified by the defect are 
very localised on the N atoms surrounding it and their closest neighbors (in the triangle 
in Fig. 2(a)), suggesting a little communication between nearby holes. A similar behavior 
is seen on the defect. Leaving an unpaired electron on each edge of the triangle, 
displays a net total magnetic moment of 3 By extension, for any triangular defect, 
the magnetic moment at ground state, Mqs^ is held by dangling bonds and should go as 
Mqs — 3 x T — 6 for T > 1. Thus, the electronic and magnetic properties of h-BN sheets 
are strongly affected by the size of the vacancy defects. 

Although the ground state of a defect is clearly M = /ig, it is interesting to 
study the impact of non-zero total magnetic moments both on its geometry and its spin 

5 



localisation. To do so, we did spin-constrained calculations by fixing the total magnetic 
moment on the whole system and minimizing the total energy, allowing full geometrical 
relaxation, in order to identify the best ordering of the local spin polarization. Figure 2(b) 
shows the relative total energy and the N-N distance at each corner (symbols defined in the 
inset) as a function of the imposed total magnetic moment — M = 0, 2, 4 and 6 /ig. As 
the total magnetic moment goes up, the total energy of the system increases monotonically 
as the homopolar N-N bonds are broken. With M = 2 the lowest-energy state with a 
non-zero magnetic moment, we observe a symmetry breaking. Two N atoms in one corner 
move away from each other, up to 2.3 A, while atoms in the two remaining pairs move 
slightly closer. Figure 2(c) shows the optimized geometry of the atoms around this defect 
with M = 2 \±b and the difference in the charge density (p^ — p±) projected on the plane of 
the monolayer, it shows that the extra spins arrange themselves predominantly in the two 
dangling N bonds arising by the opening of one corner. As the total magnetic moment is 
further increased (Fig. 2(b)) the remaining two pairs are broken in sequence, bringing the 
interatomic distance between corner-sharing N to the lattice value for second neighbor N 
atoms, i.e., 2.51 A. As the triangular defect size T 2 goes up, so does the number of unpaired 
electrons on the edges, but the N-pairing at the three corners remains. Thus, to ensure spin 
localisation the total magnetic moment that can be imposed on these triangular defects is 
M = Mqs + 2 (+4, +6), irrespective of the size. 

If a change in the total magnetic moment can affect the geometrical reconstruction, will 
a structural deformation impact the total magnetic moment? In other words, is it possible 
to modify the magnetic moment of a defect by imposing a structural deformation? Given 
that increasing M induces symmetry breaking, it is natural to first look at the application of 
an uniaxial tensile strain, either along the armchair or the zigzag direction, that would have 
the same impact. From the results on the effect of M on the reconstruction at the edge of 
the triangular defect (Fig. 2(b)), we can expect that a tensile strain in the zigzag direction 
would lead first to the rupture of the upper corner of the triangle since the N-N bond in 
this corner is oriented in the zigzag direction. Such a deformation would imply a transition 
from Mqs +2 Pb- Similarly, an uniaxial strain in the armchair direction would affect 
first the two N-N bonds in the lower corners of the triangle, with an expected transition 
from Mqs ~^ +4 /jlb- 

While the relation between an external strain and magnetic state appears relatively clear, 
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FIG. 3. (Color online), (a) Evolution of the total energy for the system with a defect as a 
function of the uniaxial deformation in the armchair direction. Open triangles, energy difference, 
AE, between magnetic states M = 4 n B and M = \ib- (b) Optimized geometry and the Wannier 
states projected on the plane of the monolayer for the system with M — 4 fiB at an uniaxial 
deformation in the armchair direction of 2.4 %. For clearness, only atoms along the borders of the 
vacancy are included, (c) Absolute and relative evolution of the N-N distance in a system with a 
defect as a function of uniaxial deformation in the armchair direction. Diamond symbols refer 
to the N-N bond in the zigzag direction, while square symbols are the mean value of both N-N 
distances as they are presented in (b). 

and it was shown, recently, that strain can indeed enhance magnetic vacancies 20 , it is the 
evolution of the total energy for the various magnetic states as a function of the deformation 
that determines whether these transformations are physically relevant or not. Figure 3(a) 
shows the evolution of the total energy as a function of the uniaxial deformation in the 
armchair direction for the system with a defect in two definite magnetic states M = \±b 
and M = 4: fiB- For clarity, we also plot the energy difference, AE 1 , between both magnetic 
states. At zero % of uniaxial deformation, as calculated before, the ground state with M = 
\±b is 1.46 eV more stable than that at M = 4 fi B . This difference reduces monotonically, 
almost linearly, as the deformation increases, showing a crossing at a deformation of about 
2.0%, at which point the M = 4 \i B becomes more stable than the one with zero magnetic 
moment (Fig. 3(a), AE line) by as much as 0.2 eV at 2.4% of deformation. 

Figure 3(b) shows the optimized geometry for a system with a total magnetic moment 
of 4 fi B at 2.4 % of deformation in the armchair direction as well as the projections of 
Wannier states on the plane of the monolayer. As can be seen, the upper bond, which is 
perpendicular to the deformation, remains closed while there is bond breaking in the bottom 
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corners. Since all other atoms maintain their original coordination, the spin magnetic charge 
density also remains localised precisely on those unpaired dangling bonds. This transition 
occurs because a large fraction of the strain is released, anharmonically, by the stretching 
of the N-N distance in the two lower corners of the triangular defect (most parallel to the 
deformation) as seen in Fig. 3(c). These extend by about 10 % at a 2.0 % of deformation 
for both states, stretching the N-N bonds in the M = \±b case to 1.8 A. The transition 
from Mqs +4 in the defect with an uniaxial strain applied along the armchair 
direction is therefore due to a very subtle balance between the overall elastic relaxation 
around the defect as the N-N bonds parallel to the deformation are broken and the cost of 
placing two spins nearby. For this mechanism, the direction of the deformation is important, 
for example in the zigzag direction the transition occurs at a greater deformation. Moreover, 
it also depends of the size of the defect: it is impossible to reproduce this phenomenon with 
uniaxial deformation along both the armchair and zigzag direction in the larger defect 
(not shown). In those cases, the cost of the elastic deformations of the whole system is higher 
than the energy gain obtained from adding pairs of spins at the corners of the triangular 
defects. However, since the arising magnetic moment remains localised, a deformation of the 
whole system around the defect is in principle not needed. The stress-induced deformations 
presented here clearly show that the N-N distance and the magnetic moment of the defect are 
strongly related to each other. It is important to point out that also a different mechanism 
which would only constrain the N-N bondings close to the defect would do the job. For 
example, this effect could also be used for molecular detection, where a molecule fixing itself 
on a defect, of various size, could induce a large enough structural deformation to induce a 
change in the spin state. This bias for N-terminated defects together with a well controlled 
irradiation time could facilitate the creation of arrays with defects, opening the door of 
a tunable spin state nanostructure. 
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